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CN i ABSTRACT 

We review some practical aspects of measuring the amplitude of variability in 'red noise' light 
curves typical of those from Active Galactic Nuclei (AGN). The quantities commonly used to 
^ [ estimate the variability amplitude in AGN light curves, such as the fractional rms variability 

■ amplitude, Fvar, and excess variance, cr'^g, are examined. Their statistical properties, relation- 



ship to the power spectrum, and uses for investigating the nature of the variability processes 
are discussed. We demonstrate that (or similarly Fvai ) shows large changes from one part 
• of the light curve to the next, even when the variability is produced by a stationary process. 

I This limits the usefulness of these estimators for quantifying differences in variability ampli- 

. tude between different sources or from epoch to epoch in one source. Some examples of the 

expected scatter in the variance are tabulated for various typical power spectral shapes, based 
on Monte Carlo simulations. The excess variance can be useful for comparing the variability 
amplitudes of light curves in different energy bands from the same observation. Monte Carlo 
O ' simulations are used to derive a description of the uncertainty in the amplitude expected be- 

tween different energy bands (due to measurement errors). Finally, these estimators are used 
^2 ' to demonstrate some variability properties of the bright Seyfert 1 galaxy Markarian 766. The 



source is found to show a strong, linear correlation between rms amplitude and flux, and to 
show significant spectral variability. 
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1 INTRODUCTION 

One of the defining characteristics of Active Galactic Nuclei 
(AGN) is that their X-ray emission is variable. X-ray light curves 
from Seyfert 1 galaxies show unpredictable and seemingly aperi- 
odic variability (Lawrence et al. 1987; M'^Hardy 1989). Such ran- 
dom variability is often referred to as noise, meaning that it is the 
result of a stochastic, as opposed to deterministic, process. In this 
context the 'noise' is intrinsic to the source and not a result of mea- 
surement errors (such as Poisson noise), i.e. the signal itself is the 
output of a noise process. 

One of the most common tools for examining AGN variability 
(and noise processes in general) is the fluctuation Power Spectral 
Density (PSD) which represents the amount of variability power 
(mean of the squared amplitude) as a function of temporal fre- 
quency (timescale^^). The high frequency PSDs of Seyferts are 
usually well-represented by power-laws over a broad range of fre- 
quencies ('P{f) oc /^°, where 'P(/) is the power at frequency /) 
with slopes a = 1 — 2 (Green, M'^Hardy & Lehto, 1993; Lawrence 
& Papadakis, 1993; Edelson & Nandra, 1999; Uttley, M'^Hardy & 
Papadakis, 2002; Vaughan, Fabian & Nandra 2003; Markowitz et 



al. 2003). Such a spectrum, with a slope a > 1 is usually called 
'red noise' (for an introduction to red noise see Press 1978). 

If Seyfert I light curves are viewed as the product of a stochas- 
tic (in this case red noise) process then the specific details of each 
individual light curve provide little physical insight. Each light 
curve is only one realisation of the underlying stochastic process, 
i.e. it is one of the ensemble of random light curves that might be 
generated by the process. Each new realisation will look different 
and these changes are simply statistical fluctuations inherent in any 
stochastic process (as opposed to changes in the nature of the pro- 
cess itself). Therefore one should expect two light curves to have 
different characteristics (such as mean and variance) even if they 
are realisations of the same process. On the other hand, data from 
deterministic processes, for example the energy spectrum of a non- 
varying source or the light curve of a strictly periodic source (such 
as a pulsar), should be repeatable within the limits set by the mea- 
surement errors. 

It is the average properties of the variability (such as the PSD) 
that often provide most insight into the driving process. For exam- 
ple, the PSD of any real red noise process cannot continue as a 
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steep power-law indefinitely to longer timescales or the integrated 
variability power would diverge. Therefore the PSDs of AGN vari- 
ability must break to a flatter index at low frequencies; the po- 
sition of such a break would represent a characteristic variability 
timescale and may yield information about the underlying driving 
process. Recent timing studies have indeed found evidence that the 
steep power-law PSDs of Seyfert Is show a flattening, or turnover, 
at low frequencies (Edelson & Nandra 1999; Uttley et al. 2002; 
Markowitz et al. 2003). 

In many cases however the data are not adequate for PSD anal- 
ysis. In these situations the variability is usually described in terms 
of the statistical moments (e.g. the sample mean and variance, etc.). 
However, due to the stochastic nature of red noise variability there 
is a large degree of randomness associated with these quantities. In 
practice this means that it is difficult to assign meaningful errors to 
the variance. This in turn makes it difficult to quantitatively com- 
pare variances, say from repeated observations of the same source 
(and thereby test whether the variability is stationary). Such an 
analysis might be desirable; it could in principle reveal changes 
in the 'state' of the source if its variability properties were found 
to evolve with time. This paper discusses this and related problems 
that are encountered when examining the variability properties of 
AGN. Particular emphasis is placed on the mathematical properties 
and implications of the inherent randomness in the variability. The 
mathematical details are well understood from the general theory 
of stochastic processes (e.g. Priestley 1981 for spectral analysis) 
but some of the practical consequences for AGN observations have 
not been discussed in detail. On the basis of simulated data some 
recipes are developed that may serve as a useful guide for observers 
wishing make quantitative use of their variability analysis without 
the recourse to e.g. extensive Monte Carlo simulations. 

The paper is organised as follows. Section|2|defines the esti- 
mators to be discussed (namely the periodogram and the variance). 
Simulated data are used to illustrate various aspects of these es- 
timators; section |3] describes methods for producing artificial red 
noise time series. Section |4] discusses the stationarity of time se- 
ries. Sections|5]and|6|discuss two sources of uncertainty associated 
with measuring variability amplitudes, the first due to the stochas- 
tic nature of the variability and the second due to flux measurement 
errors. Section0gives an example using a real XMM-Newlon ob- 
servation of Mrk 766. Finally, a brief discussion of these results is 
given in section|8|and the conclusions are summarised in section|9| 



2 ESTIMATING THE POWER SPECTRAL DENSITY 

The PSD defines the amount of variability 'power' as a function of 
temporal frequency. It is estimated by calculating the periodogram^ 
(Priestley 1981; Bloomfield 2000). 

For an evenly sampled light curve (with a sampling period 
AT) the periodogram is the modulus-squared of the Discrete 
Fourier Transform (DFT) of the data (Press et al. 1996). For a light 
curve comprising a series of fluxes Xi measured at discrete times 
t,:(i = l,2,...,7V): 
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(1) 



at N/2 evenly spaced frequencies fj = j/NAT (where j = 
1,2,..., 7V/2), fN/2 = l/2Ar is the Nyquist frequency, /Nyq. 
Note that it is customary to subtract the mean flux from the 
light curve before calculating the DFT. This eliminates the zero- 
frequency power. The periodogram, P{fj), is then calculated by 
choosing an appropriate normalisation A (see AppendixIXIfor more 
on periodogram normalisations). For example 
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If the time series is a photon counting signal such as normally 
encountered in X-ray astronomy, and is binned into intervals of 
AT, the effect of Poisson noise is to add an approximately con- 
stant amount of power to the periodogram at all frequencies. With 
the above normalisation this constant Poisson noise level is 2x (as- 
suming the light curve is not background subtracted). 



2.1 Statistical properties of the periodogram 

The periodogram of a noise process, if measured from a single time 
series, shows a great deal of scatter around the underlying PSD. In 
particular, the periodogram at a given frequency [P(/)] is scat- 
tered around the PSD ['P{f)] following a distribution with two 
degrees of freedom (van der Klis 1989): 



P{f) = P{f)xl/2, 



(3) 



where X2 is a random variable distributed as with two degrees of 
freedom, i.e. an exponential distribution with a mean and variance 
of two and four, respectively. The periodogram is distributed in this 
way because the real and imaginary parts of the DFT are normally 
distributed for a stochastic process^ (section 6.2 of Priestley, 1981; 
Jenkins & Watts 1968). The expectation value of the periodogram 
is equal to the PSD but its standard deviation is 100 per cent, lead- 
ing to the larger scatter in the periodogram (see Fig.0. See Leahy 
et al. (1983), van der Khs (1989), Papadakis & Lawrence (1993), 
Timmer & Konig (1995) and Stella et al. (1997) for further discus- 
sion of this point. 

When applied to real data the periodogram is an inconsistent 
estimator of the PSD, meaning that the scatter in the periodogram 
does not decrease as the number of data points in the light curve 
increases (Jenkins & Watts 1968). In order to reduce this scatter 
the periodogram must be smoothed (averaged) in some fashion. As 
the number of data points per bin increases (either by binning over 
frequencies or averaging over many data segments) the scatter in 
the binned periodogram decreases, i.e. the averaged periodogram 
is a consistent estimator of the PSD (see Papadakis & Lawrence 
1993 and van der Klis 1997 for more on binned periodogram esti- 
mates). A further point is that periodograms measured from finite 
data tend to be biased by windowing effects which further compli- 
cate their interpretation (van der Klis 1989; Papadakis & Lawrence 
1993; Uttley et al. 2002 and see below). 



^ Following Priestley (1981) the term "periodogram" is used for the dis- 
crete function P{fj), which is an estimator of the continuous PSD P{f). 
The periodogram is therefore specific to each realisation of the process, 
whereas the PSD is representative of the true, underlying process. 



^ The DFT at the Nyquist frequency is always real when A'^ is even so the 
periodogram at this frequency is distributed as xf, i e. with one degree of 
freedom. 
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Time Frequency 

Figure 1. Simulated time series (left) and their periodograms (right). The 
upper panel shows a 'flicker noise' time series which has a PSD. The 
lower panel shows a 'random walk' time series with a. PSD. Note the 
large scatter in the periodogram (dots) around the underlying PSD (solid 
line). It is clear that the time series with the steeper PSD shows more power 
in long-term variabiUty while the time series with the flatter PSD shows 
relatively more power in short term variability (flickering). The two series 
were generated using the same random number sequence. 



2.2 Integrated power 

The integral of the PSD between two frequencies (/i and /2) yields 
the contribution to the expectation value of the ('true') variance 
due to variations between the corresponding timescales (I//1 and 
I//2 ). This result follows from Parseval's theorem (see e.g. van der 
Khs 1989) 

(s') = f' Vim. (4) 
-I fi 

Correspondingly, for a discrete time series the integrated peri- 
odogram yields the observed variance for that particular realisation 

JV/2 

5'=5^P(/.)A/, (5) 

where A/ is the frequency resolution of the DFT (A/ = 
1/NAT). The total variance of a real light curve is equal to its 
periodogram integrated over the frequency range fi — 1 /NAT to 

/Nyq = l/2Ar. 

The sample variance (which will differ from observation to 
observation) is given by: 

JV 

^' = iV^E(-'-^)'' (6) 

i = l 

where x is the arithmetic mean of Xi. In the limit of large TV these 
two variance estimates are identical. The normalised variance"^ is 
simply S^/x^. 



In AGN studies normalised quantities are often used in preference to ab- 
solute quantities as they are independent of the flux of a specific source. 
This means that, in principle, normalised amplitudes can be used to com- 
pare sources with different fluxes. 



3 SIMULATING RED NOISE LIGHT CURVES 
3.1 Algorithms 

In order to elucidate the properties of the variance of red noise data, 
random light curves were generated from power-law PSDs similar 
to those of AGN. Fig. Q shows two artificial time series and their 
periodograms. It is worth reiterating that the large scatter in the 
periodograms is an intrinsic property of stochastic processes - it 
does not depend on the number of data points and is not related to 
Poisson noise in the data. 

These artificial time series were produced using the algorithm 
of Timmer & Konig (1995). This generates random time series with 
arbitrary broad-band PSD, correctly accounting for the intrinsic 
scatter in the powers (i.e. equation|3}- Other methods of generating 
random light curves include the related 'summing of sines' method 
(Done et al. 1992). Note that it is not correct to randomise only the 
phases of the component sine functions, their amplitudes must also 
be randomised. Otherwise this method does not account for this 
intrinsic scatter in the powers. Shot-noise models can produce red 
noise time series with certain PSD shapes (see Lehto 1989). There 
also exist various mathematical tricks for producing data with spe- 
cific power-law PSD slopes. Data with a a = 1 PSD (often called 
'flicker noise') can be generated using the half-integral method out- 
lined in Press (1978), while a = 2 ('random walk') data can be 
generated using a first-order autoregressive process (AR[1]), essen- 
tially a running sum of Gaussian deviates (see Deeming 1970 and 
Scargle 1981 for more on such methods). The method of Timmer 
& Konig (1995) is used below as this can generate time series from 
an arbitrary PSD and is computationally efficient. 



3.2 Simulating 'realistic' data 

Some caution should be applied when using these routines to pro- 
duce artificial time series. As mentioned briefly in the previous sec- 
tion, periodograms measured from real data tend to be biased by 
windowing effects. For uninterrupted but finite observations data 
of red noise processes the most important of these effects is 'red 
noise leak' - the transfer of power from low to high frequencies 
by the lobes of the window function (see e.g. Deeter & Boynton 
1982; van der Klis 1997). If there is significant power at frequen- 
cies below the lowest frequency probed by the periodogram (i.e. on 
timescales longer than the length of the observation) this can give 
rise to slow rising or falling trends across the light curve. These 
trends contribute to the variance of the light curve. Thus variability 
power 'leaks' into the data from frequencies below the observed 
frequency band-pass. The degree to which this occurs, and the re- 
sultant bias on the measured periodogram, depend on the shape of 
the underlying PSD and the length of the observation (Papadakis & 
Lawrence 1995; Uttley et al. 2002). For flat PSD slopes (a < 1.5) 
the amount of leakage from low frequencies is usually negligible. 

Since AGN light curves usually contain significant power on 
timescales longer than those probed (see section l4m the effects 
of red noise leak must be included in simulations of AGN light 
curves. This can be achieved by using the Timmer & Konig (1995) 
algorithm to generate a long light curve from a PSD that extends to 
very low frequencies and then using a light curve segment of the re- 
quired length. Data simulated in such a fashion will include power 
on timescales much longer than covered in the short segment. The 
effects of measurement errors (e.g. Poisson noise) can be included 
in the simulations using standard techniques (e.g. Press et al. 1996). 
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4 STATIONARITY 

A stationary process is one for which the statistical properties (such 
as mean, variance, etc.) do not depend on time. Fig.|2|shows an ar- 
tificial red noise time series together with its mean and variance 
measured every 20 data points. The simulation was produced from 
a single, well-defined PSD (which did not vary). The process is 
therefore stationary. It would have been reasonable to expect the 
resulting time series (a realisation of the process) to appear sta- 
tionary. This is not the case however; both the mean and variance 
change with time (panels 2 and 3). This has nothing whatsoever to 
do with measurement errors - the simulation has zero errors. This 
simulation demonstrates that, when dealing with red noise, fluctu- 
ations in variance are not sufficient to claim the variability process 
is non-stationary. 

As the purpose of time series analysis is to gain insight into 
the process, not the details of any specific realisation, a more robust 
approach is needed to determine whether these data were produced 
by a time-stationary process or a non- stationary process. It would 
be more insightful to consider whether the expectation values of 
the characteristics (such as the variance) are time- variable. The ex- 
pectation values should be representative of the properties of the 
underlying process, not just any one realisation. See section 1.3 of 
Bendat & Piersol (1986) for a discussion of this point. 



4.1 Weak non-stationarity 

For a process with a steep red noise PSD (a ^ 1), the integrated 
periodogram will diverge as / — > 0. This means that (following 
equation|4} the variance of a red noise time series with a steep PSD 
will diverge with time. In this case there is no well-defined mean; 
Press & Rybicki (1997) describe this form of variability as 'weakly 
non-stationary.' For time series with power spectra flatter than this 
the variance converges as / — > 0, thus for a white noise process 
with a flat PSD (a — 0), the variance will converge as the obser- 
vation length increases, and there will be a well-defined mean on 
long timescales. 

Of course, for any real process the PSD must eventually flat- 
ten such that the power does not diverge (i.e. a < 1 on sufficiently 
long timescales). Thus, weak non-stationarity is entirely due to ob- 
servations sampling only the steep part of the PSD of a source. 
But in AGN this flattening occurs on timescales much longer than 
those probed by typical observations. For instance, XMM-Newton 
observations of AGN typically last for ~ few x 10* s whereas in 
many objects the PSD is steep until > IQp s and in some cases 
probably much longer (Edelson & Nandra 1999; Uttley et al. 2002; 
Markowitz et al. 2003). Therefore on the timescales relevant for 
most X-ray observations, AGN light curves should be considered 
weakly non-stationary. 



4.2 Stochasticity 

Fluctuations in the statistical moments (such as mean and variance) 
of a light curve are intrinsic to red noise processes. Therefore, even 
in the absence of measurement errors (e.g. no Poisson noise) the 
means and variances of two light curves produced by exactly the 
same process can be significantly different. This can be seen in 
Fig. |2| (panels 2 and 3), where each 20 point segment of the light 
curve shows a different mean and variance. These random fluctu- 
ations in variance are however governed by the normal statistical 
rules of noise processes and can thus be understood in a statistical 
sense. 
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Figure 2. Panel 1: Simulated red noise time series (with a /^^ PSD) with 
N = 2800 points. Panel 2 and 3: mean and variance measured from seg- 
ments of 20 points (calculated using equation 1^. The variances follow a 
distribution of the form shown in the bottom panel of Fig.|3] (Note the vaii- 
ance is plotted on a logarithmic scale.) Panel 4: averaged vaiiance measured 
by binning the individual variances into groups of 20 consecutive estimates. 
The errors are the standard error on the mean (equation 4.14 of Bevington 
& Robinson 1992). These averaged variances are consistent with staying 
constant. In other words, although the instantaneous value of the variance 
fluctuates, its expectation value is consistent with constant (i.e. a station- 
ary process). Panel 5: fractional rms amplitude (yj S"^ /x^) measured from 
segments of 20 points. Panel 6: averaged fractional rms amplitude measured 
by binning the individual amplitudes into groups of 20. The fractional am- 
plitude is anti-correlated with the light curve because {S^) is constant but 
-Fvar is normalised by the light curve flux. 



Any given series is only one realisation of the processes and 
its periodogram will show the scatter predicted by equation|3| The 
integrated periodogram (which gives the variance; equation|5} will 
therefore be randomly scattered around the true value for the PSD 
of the process. The variance in a specific time series is given by 



N/2 



NAT 



(7) 



i.e. the variance of a given realisation is a sum of X2 distributions 
weighted by the PSD''. (This assumes biases such as red noise 
leak are not significant. If this is not true then these biases will 

* As noted earlier, the periodogram at the Nyquist frequency is actually 
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Figure 3. Distribution of variances in time series witli tliree different PSD 
shapes: f'^ (top), (middle) and (bottom). Each distribution is de- 
rived from 50,000 realisations. As the PSD gets steeper the distribution of 
variances becomes less Gaussian and more like a distribution with a low 
effective degrees of freedom. 

further distort the distribution of variances.) It is the expectation 
value of the variance that is representative of the integrated power 
in the PSD, and thus the average amplitude of the variability pro- 
cess (eqn.|3). Thus, while the expectation value of the variance is 
equal to the integrated PSD, each realisation (time series) of the 
same process will show a different variance even if the parent vari- 
ability process is stationary. This is particularly important for steep 
PSD time series (i.e. weakly non-stationary data) since the vari- 
ance does not converge as more data are collected. Only if the time 
series spans timescales on which the integrated power converges 
at low frequencies (i.e. q < 1) will the variance converge as the 
length of time series increases. 

distributed as Xi for sven N. However, for large this will make a negli- 
gible difference to the sum (cf. equation 2.9 of van der Klis, 1989). 



These points are illustrated by Fig.|3| which shows the distri- 
butions of variances in random time series with three different PSD 
slopes (a — 0, 1, 2). This plot was produced by generating 50,000 
random time series (each 100 points long) for each PSD slope and 
measuring the variance of each one. The scatter in the variance is 
entirely due to random fluctuations between different realisations 
because the PSD normalisation was kept fixed and no instrumental 
noise was added. The shape of the distribution of variances can be 
seen to depend on the PSD slope. 

Consider a white noise process (a — 0, i.e. 7-'(/) = const). 
The periodogram of each realisation is randomly scattered around 
its flat PSD. Following equation |7| the variance is simply the sum 
of the N/2 xi -distributed powers in the periodogram, and these 
are evenly weighted (PSD is constant). The sum of N/2 xi distri- 
butions follows a Xat distribution. This tends to a normal distribu- 
tion as A'^ increases. Thus the variance of a white noise process is 
approximately normally distributed (as can be seen in Fig.|3j and 
converges as A'^ increases. The fractional standard deviation of %% 
is given by ^2/N , so for the 100 point light curves used the (Ict) 
fractional width of the variance distribution is ~ 14.1 per cent, in 
agreement with the simulations (Fig.|3] top panel). 

For time series with steeper PSDs the lower frequency peri- 
odogram points contribute more strongly to the sum than the higher 
frequency points. The variance of such a time series is therefore 
dominated by a few low frequency powers^ and thus resembles a 
X^ distribution with low 'effective degrees of freedom.' The distri- 
bution of variances in red noise data is dependent on the underlying 
PSD and is, in general, non-Gaussian (Fig.|3j- The fractional stan- 
dard deviation of xi5 is \/2/v, which tends to unity as the PSD 
gets steeper (i.e. as effective v 2). Thus the largest fluctuations 
in variance (up to a limit of ~ 100 per cent rms) are expected to 
result from very steep PSD slopes. 



5 INTRINSIC SCATTER IN VARIANCE 

As discussed above, when examining AGN light curves one should 
expect random changes in the mean and variance with time (be- 
tween segments of a long observation or between observations 
taken at different epochs). This is true even if the measurement 
errors are zero and does not depend on the number of data points 
used (due to the weak non-stationarity). However, it is also possi- 
ble that the underlying process responsible for the variability itself 
changes with time (e.g. the PSD changes), in which case the vari- 
ability is non-stationary in a more meaningful sense - 'strongly 
non-stationary.' Such changes in the variability process could pro- 
vide insight into the changing physical conditions in the nuclear 
regions. On the other hand the random changes expected for a red 
noise process yield no such physical insight. The question thus 
arises: how does one tell, from a set of time series of the same 
source, whether they were produced by a strongly non-stationary 
process? In other words, is it possible to differentiate between dif- 
ferences in variance caused by real changes in the variability pro- 
cess (physical changes in the system), and random fluctuations ex- 
pected from red noise (random nature of the process)? 

If the process responsible for the variability observed in a 



^ The windowing effects mentioned above mean that fluctuations at powers 
above and below the frequency range of the periodogram may also affect the 
variance. 
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given source is stationary tlien its PSD is constant in time. Tlie ex- 
pectation value of the absolute (un-normalised) variance will there- 
fore be the same from epoch to epoch, but the individual variance 
estimates will fluctuate as discussed in section l4!2l This makes it 
difficult to judge, from just the variances of two light curves taken 
at two epochs, whether they were produced by a stationary process. 
Given sufficient data it is, however, possible to test whether the ex- 
pectation values of the variance (estimated from an ensemble of 
light curves) at two epochs are consistent with a stationary process. 

5.1 Comparing PSDs 

The methods most frequently employed involve comparing the 
PSDs (estimated from the binned periodogram) at different epochs. 
If the PSDs show significant differences (at a given confidence 
level) the variability process can be said to be strongly non- 
stationary. As an example of this, the PSDs of X-ray binaries 
evolve with time, and the way in which the variability properties 
evolve provides a great deal of information on the detailed work- 
ings of these systems (see e.g. Belloni & Hasinger 1990; Uttley & 
M'^Hardy 2001; Belloni, Psaltis & van der Klis 2002; Pottschmidt 
et al. 2002). 

Papadakis & Lawrence (1995) suggested a method suitable for 
testing whether large AGN datasets display evidence for strongly 
non-stationary variability. Again this method works by comparing 
the PSDs from different time intervals, in this case by determining 
whether the differences between two periodograms are consistent 
with the scatter expected based on equation|3| In particular, they de- 
fine a statistic s based on the ratio of two normalised periodograms. 
If s deviates significantly from its expected value for stationary data 
(if (s) =0) then the hypothesis that the data are stationary can be 
rejected (at some confidence level). 

5.2 Comparing variances 

A different approach is compare variances Si derived from Al ob- 
servations of the same source (either segments of one long obser- 
vation or separate short observations). In order to test whether the 
Si differ significantly (i.e. more than expected for a red noise pro- 
cess) a measure of the expected scatter is required. This error could 
be obtained directly from the data (by measuring the standard de- 
viation of multiple estimates) or through simulations (based on an 
assumed PSD shape)^. 

5.2.1 Empirical error on variance 

An empirical estimate of the mean and standard deviation of the 
variance can be made given M non-overlapping data segments. The 
M segments each yield an estimate of the variance^. Si. Each of 
these is an independent variable of (in general) unknown but identi- 
cal distribution (unless the process is strongly non-stationary). The 
central limit theorem dictates that the sum of these will become 

^ It is assumed that the data segments being compared have identical sam- 
pling (same bin size and observation length). Every effort should be made 
to ensure this is the case, e.g. by clipping the segments to the same length. 
The variances will then be calculated over the same range of timescale (fre- 
quencies). As the variance can increase rapidly with timescale in red noise 
data this is most important for steep PSD data such as AGN light curves. 
^ Ideally each segment should contain at least A'^ > 20 data points in order 
to yield a meaningful variance. 
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Figure 4. The average rms amplitude (a = V 5^) as a function of flux for 
the simulated light curve shown in Fig.|2] The individual rms estimates were 
sorted by flux and binned to M = 20 estimates per bin. Errors correspond 
to the error on the mean value. The amplitude is constant with flux. 

normally distributed as M increases. Therefore by averaging the 
Al variance estimates it is possible to produce an averaged vari- 
ance ((S^)) and assign an error bar in the usual fashion (e.g. equa- 
tion 4.14 of Bevington & Robinson 1992). This gives a robust esti- 
mate of the variance and the standard deviation of the Al variances 
around the mean gives an estimate of the uncertainty on the mean 
variance. 

If several sets of data segments are acquired it is therefore pos- 
sible to compare the mean variance of each set statistically (since 
each has an associated uncertainty). For example, with two long 
XMM-Newton observations of the same source, taken a year apart, 
one could measure the variance for each observation (by breaking 
each into short segments and taking the mean variance of the seg- 
ments). Thus it would be possible to test whether the variability 
was stationary between the two observations. This method of esti- 
mating the mean and standard deviation of the variance requires a 
large amount of data. Of order iV x M = 20 x 20 = 400 data 
points are needed to produce a single well-determined estimate of 
the mean variance and its error. A typical XMM-Newton observa- 
tion of a bright Seyfert 1 galaxy (~ 40 ks duration) is only likely 
to yield enough data for one estimate of the mean variance. Thus 
this method is suitable for testing whether the mean variance has 
changed from observation to observation. 

Fig. |2| (panel 4) demonstrates this empirically derived mean 
variance and its error bar on a long, simulated time series. These 
data were produced by calculating the variances 5*^ in bins of 
TV = 20 data points (panel 3) and then averaging Al = 20 vari- 
ances to produce a mean variance with error bar (panel 4). These 
averaged variances are consistent with constant, as expected; fit- 
ting these data with a constant gave xi ~ 0.84. Figure |4] shows 
the rms amplitude is constant with flux. These tests indicate that 
the integrated PSD is consistent with being constant with time; the 
variance does not change significantly from epoch to epoch (or as 
function of flux), as expected for a stationary process. 

5.2.2 Estimating the error on the variance through simulations 

The advantage of the above method is that it requires no assump- 
tion about the shape of the PSD, The drawback is that it requires a 
substantial amount of data to produce a single, robust variance esti- 
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mate. An alternative approach is to estimate the standard deviation 
of the variances Si based on simulations. 

Given an assumed shape for the PSD it is possible to calcu- 
late the distribution of variances expected for a stationary process 
(see section l4!2l . Some example distributions are shown in Fig.|3| 
which clearly demonstrates how the distribution depends on the 
slope of the PSD. The distribution becomes more normal at flat- 
ter slopes and more asymmetric at steep slopes. For a given PSD 
shape these distributions are well defined (by eqn. and can be 
computed through Monte Carlo simulations. This makes it possible 
to estimate limits within which one would expect the variance to be 
distributed if the process is stationary. 

The two primary factors that affect the distribution of variance 
are the PSD shape and the sampling of the light curve (the length of 
the data segments in the case of contiguously binned light curves). 
TableQ gives the expected confidence limits for four different PSD 
shapes and five different lengths for the data segments. These val- 
ues were computed by simulating one very long light curve with 
the assumed PSD shape and breaking it into 1000 separate seg- 
ments (of specified length). The variance within each segment was 
measured and the distribution of the 1000 variances was calculated. 
The 90 per cent confidence interval was calculated by finding the 
5th and 95th percentiles of the variance distribution (in general 
these upper and lower bounds will differ because the distribution 
is asymmetric). The numbers given in the table are the boundaries 
of the 90 and 99 per cent confidence regions estimated by aver- 
aging the results from 50 runs. The limits are given in terms of 
±A log(S^) because they are multiplicative. That is, from a partic- 
ular realisation the variance is expected to be scattered within some 
factor of the true variance (for which the absolute normalisation 
is irrelevant). The factors are tabulated in terms of their logarithms 
(since multiplicative factors in linear-space become additive offsets 
in log-space). 

The PSD used for the simulations was chosen to match that ex- 
pected for AGN, i.e. a steep power-law at high frequencies (with a 
slope of a = 1.0, 1.5, 2.0, 2.5) breaking to a flatter slope (q = 1.0) 
at low frequencies. The frequency of the break was fixed to be 
10^'^, in other words the break timescale was 1000 times the bin 
size. (The absolute size of the time bins is arbitrary in the simula- 
tions. When comparing the simulated results to real data sampled 
with e.g. 25 s time resolution, the break timescale in the simulated 
PSD is thus 25 ks.) 

The numbers given in the table provide an approximate pre- 
scription for the expected scatter in the variance of a stationary 
process with a red noise PSD similar to that of AGN. The simu- 
lated light curve shown in Fig.|2|was used to demonstrate the use 
of this table. In this case the PSD is know to have a slope a = 2, 
and the variances (shown in panel 3 of Fig.|2} were calculated ev- 
ery 20 points. Therefore, the 90 interval for the expected variance 
is given by log(S'^)to'7i- Taking the mean variance as the expec- 
tation value for S^ this translates to = 59.9 (11.7 - 168.8). 
The interval boundaries were calculated by converting the logarith- 
mic value into a linear factor and multiplying by the sample mean 
(assumed to represent the true variance). This interval is shown on 
Fig.|5|by the dotted lines. The corresponding 99 per cent confidence 
interval is also marked. 

As expected the individual variances fall within the expected 
region. However, the 90 per cent region spans an order of magni- 
tude in variance. Thus even order of magnitude differences in vari- 
ance between short sections of a light curve are to be expected and 
do not necessarily indicate that the underlying process is not sta- 
tionary. Subtle changes in the PSD will thus be difficult to detect 



Table 1. Expected scatter in variance estimates. The 90 and 99 per cent 
intervals are presented in terms of it A log(5^). (The 99 per cent interval 
is given in bold.) The boundaries were calculated from Monte Carlo simu- 
lations of light curves. The PSD was chosen to be a broken power-law with 
a slope of a = 1 below the break (at a frequency 10^'^) and a slope above 
the break ofa = 1.0, 1.5, 2.0, 2.5. The simulated data segments were cho- 
sen to be 10, 20, 50, 100, 1000 points long (the timescale of the break in 
the PSD being at 1000, in arbitrary units). 
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Figure 5. Variance of the simulated data shown in Fig.|2|(panel 3) with the 
90 (dotted line) and 99 per cent (dashed line) confidence intervals marked 
(as calculated in section l7.2.2r Clearly the variances fall within these limits, 
as expected for a stationary process. The solid line marks the mean variance. 



by examining the raw variances as the intrinsic scatter is so large. 
Such changes could be revealed by comparing averaged variances 
or comparing the PSDs as described above. 



6 EFFECT OF MEASUREMENT ERRORS 
6.1 Excess variance and Fvar 

The datasets considered thus far have been ideal, in the sense that 
they are free from flux uncertainties. In real life, however, a light 
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curve Xi will have finite uncertainties CTcrr.i due to measurement er- 
rors (such as Poisson noise in the case of an X-ray photon counting 
signal). These uncertainties on the individual flux measurements 
will contribute an additional variance. This leads to the use of the 
'excess variance' (Nandra et al. 1997; Edelson et al. 2002) as an 
estimator of the intrinsic source variance. This is the variance after 
subtracting the contribution expected from measurement errors 



2 o2 



where crfrr is the mean square error 
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The normalised excess variance is given by cr^xs — ""xs/^^ ™d 
the fractional root mean square (rms) variability amplitude (Fvarl 
Edelson, Pike & Krolik 1990; Rodriguez-Pascual et al. 1997) is the 
square root of this, i.e. 



S2 



(10) 



The statistic JVar is often chosen in preference to (Jj^xs ^1" 
though the two convey exactly the same information. Fvar is a 
linear statistics and can therefore give the rms variability ampli- 
tude in per centage terms. The choice of whether to quote Fvar 
or CTnxs is usually purely one of presentation. It is worth not- 
ing that the Monte Carlo results given in section 15.2.21 to esti- 
mate the expected scatter on the variance, can also be applied 
to its square root. The expected boundaries of the confidence re- 
gion of the logarithm of the rms is approximately half those of 
the variance. Specifically, Alog((7) ~ Alog(S'^)/2 and similarly 

Alog(Fvar) ~ Al0g(0-|,xs)/2. 



6.2 Spectral variability 

An X-ray light curve of an AGN can be split into different energy 
bands. The light curves in each band will be strictly simultaneous 
and can be used to test whether the X-ray variability is a function 
of energy. For example, one might examine the ratio of a soft band 
light curve to a hard band light curve. The statistical significance 
of any variations in the ratio can be quantified by propagating the 
measurement errors and applying an appropriate test, such as the 
test of the constant ratio hypothesis*. If the ratio shows varia- 
tions greater than those expected from the errors then the two light 
curves are intrinsically different and the source does indeed show 
spectral variability. Such changes in the energy spectrum with time 
can in principle provide valuable clues to the nature of the X-ray 
source. This test does not provide any quantitative description of 
the spectral variability. 

Another tool for investigating spectral variability is the rms 
spectrum, i.e. the rms variability amplitude (or Fvar) as a func- 
tion of energy. See e.g. Inoue & Matsumoto (2001), Edelson et al. 
(2002), Fabian et al. (2002) and Schurch & Warwick (2002) for 
some examples of rms spectra from AGN. However, when exam- 
ining rms spectra it is often not clear whether changes in the am- 
plitude with energy reflect real energy-dependence of the intrinsic 
variability amplitude or are caused by random errors in the light 
curves. The finite measurement errors on the individual fluxes (e.g. 



° This does, of course, assume the light curves have been binned suffi- 
ciently for the error bars to be approximately Gaussian. 



due to Poisson noise) will introduce some uncertainty in the esti- 
mated rms amplitudes. An estimate of this uncertainty would help 
answer the question posed above, namely whether features in rms 
spectra are the result of random errors in the data or represent spec- 
tral variations intrinsic to the source. 

The problem of how to assess the uncertainty on the excess 
variance (or Fvar) is a long-standing one (e.g. Nandra et al. 1997; 
Turner et al. 1999; Edelson et al. 2002). The standard error formu- 
lae presented in the literature (e.g. Turner et al. 1999; Edelson et al. 
2002) are formally valid in the case of un-correlated Gaussian pro- 
cesses. Typically AGN light curves at different X-ray energies are 
strongly correlated and are not Gaussian. However, when searching 
for subtle differences in amplitude between simultaneous and cor- 
related light curves it may be more useful to have an indication of 
the uncertainty resulting from the finite flux errors. 



6.2.1 Uncertainty on excess variance due to measurement errors 

A Monte Carlo approach was used to develop a prescription of the 
effect of measurement errors on estimates of Fvar (and cr|[xs). A 
short red noise light curve was generated. Poisson noise was added 
(i.e. the individual flux measurements were randomised following 
the Poisson distribution) and the excess variance was recorded. The 
fluxes of the original light curve were randomised again and the ex- 
cess variance recorded, this was repeated many times. The distribu- 
tion of excess variances was then used to determine the uncertainty 
in the variance estimate caused by Poisson noise. Full details of the 
procedure are given in Appendix IbI 

For these simulations it was found that the error on cr^xs de- 
creases as the S/N in the light curve is increased according to: 



err(o-Nxs) 




See appendix B for details of this equation and its equivalent in 

terms of Fvar. 

As this only accounts for the effect of flux measurement er- 
rors (such as Poisson noise) in a given light curve it can be used 
to test whether two simultaneously observed light curves of the 
same source, but in different bands, show consistent amplitudes. 
A demonstration of this using real data is given in the following 
section. This uncertainty does not account for the random scatter 
intrinsic to the red noise process, therefore the absolute value of 
the rms spectrum will change between realisations (i.e. from epoch 
to epoch). But if a source shows achromatic variability then the 
values of Fvar calculated in each energy band (at a given epoch) 
should match to within the limits set by the Poisson noise (i.e. the 
fractional rms spectrum should be constant to within the uncertain- 
ties given by the above equation). Differences in Fvar significantly 
larger than these would indicate that the source variability ampli- 
tude is a function of energy. This would then mean the PSD ampli- 
tude/shape is different in different energy bands, or there are multi- 
ple spectral components that vary independently. 

The above uncertainty estimates can be used to test the hy- 
pothesis that the source variability is achromatic. If significant dif- 
ferences between energy bands are detected (as in the case of Mrk 
766 presented below) then these errors should not be used to fit 
the rms spectrum. The assumption that the differences are due only 
to measurement errors is no longer the case. In such situations the 
light curves in adjacent energy bands are likely to be partially cor- 
related and so -fitting of the rms spectrum is not appropriate. The 
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Figure 6. Top panel: 0.2-10.0 keV pn light curve of Mrk 766 (with 25 s 
bins, in units of ct s~^). Panel 2 and 3: mean count rate and excess vari- 
ance measured from segments of 20 points. Panel 4: averaged normalised 
excess variance measured by binning the individual variance estimates into 
groups of 20. This average variance is inconsistent with constant. Panel 5: 
fractional rms amplitude measured from segments of 20 points. Panel 6: 
averaged fractional rms amplitude measured by binning the individual am- 
plitudes into groups of 20. This average fractional amplitude is consistent 
with constant. This contrasts with the situation shown in Fig.l2l 



2x10* 



4x10* 6x10* 
Time (s) 



BxlO* 



lo-" 



Figure 7. Excess variance of the Mrk 766 data shown in Fig.|^(panel 3) 
with the 90 (dotted line) and 99 per cent (dashed line) confidence inter- 
vals marked (as calculated in section l5.2.2l . The variances fall within these 
limits, as expected for a stationary process. The soUd line marks the mean 
variance. 



7.1 Observation details 

Mrk 766 was observed by XMM-Newton (Jansen et al. 2001) over 
the period 2001 May 20 - 2001 May 21 (rev. 265). The present 
analysis is restricted to the data from the pn European Photon Imag- 
ing Camera (EPIC), which was operated in small window mode. 
Extraction of science products from the Observation Data Files 
(ODFs) followed standard procedures using the XMM-Newton Sci- 
ence Analysis System (SAS) v5.3.3. Source data were extracted 
from a circular region of radius 35 arcsec from the processed image 
and only events corresponding to patterns 0-4 (single and double 
pixel events) were used. Background events were extracted from 
regions in the small window least effected by source photons, these 
showed that the background rate increased dramatically during the 
final ~ 1.5 X 10^ s of the observation. This section of the data 
was excluded, leaving 1.05 x 10"" s of uninterrupted data. The light 
curves were corrected for telemetry drop outs (less than 1 per cent 
of the total time) and background subtracted. The errors on the light 
curves were calculated by propagating the Poisson noise. 



differences in excess variance will be a combination of intrinsic dif- 
ferences and measurement errors. Their uncertainty will therefore 
be more difficult to quantify. 



7 CASE STUDY: AN XMM-Newton OBSERVATION OF 
Mrk 766 

In this section a long (~ 10^ s) XMM-Newton observation of the 
bright, variable Seyfert I galaxy Markarian 766 is used to illus- 
trate the points discussed above. The data were obtained from the 
XMM-Newton Data Archive^. Details of the observation are dis- 
cussed in Mason et al. (2003) and an analysis of the PSD is pre- 
sented in Vaughan & Fabian (2003). 



7.2 Stationarity of the data 

The broad band (0.2-10 keV) light curve extracted from the pn is 
shown in Fig. |6| (panel 1). As was the case for the simulated data 
shown in Fig|2| the mean and variance (calculated every 20 data 
points) show changes during the length of the observation (pan- 
els 2 and 3). The expected range for the excess variance, calculated 
using the results of section 15.2.21 (and assuming a PSD slope of 
a = 2.0), is marked in figure Q Fig. |8| shows the same data in 
terms of normalised excess variances. Neither of these show fluc- 
tuations larger than expected for a stationary process. But given the 
large expected scatter this is a rather insensitive test. In the case of 
the Mrk 766 light curve however, there are sufficient data to exam- 
ine variations of the average variance with time, allowing a more 
sensitive test for non-stationarity. 



" These correspond to 'instantaneous' estimates of the source variance on 
' |http : //xmm . vilspa ■ esa ■ es] timescales of 50 — 500 s. 



10 Vaughan et al. 



5 10 15 20 25 



Time (s) 

Figure 8. As for Fig.lTlbut using the normalised excess variance of the Mrlc 
766 data. 



By averaging the excess variance estimates (in time bins con- 
taining 20 excess variance estimates) significant clianges in the 
variance with time are revealed (panel 4). This contrasts with the 
simulated data shown in Fig.|2| The binned excess variance is in- 
consistent with a constant hypothesis: fitting with a constant gave 
= 23.1 for 9 degrees of freedom (dof), rejected at 99 per cent 
confidence. The average variance is therefore changing with time, 
indicating the variability is strongly non-stationary. 

A careful inspection of Fig.|6|(panels 3 and 4) shows the indi- 
vidual variance estimates have a tendency to track the source count 
rate. This is difficult to discern from the individual variances (panel 
3), due to the larger intrinsic scatter, but much clearer in the av- 
eraged variances (panel 4). This can be seen clearly in Fig.|9|(top 
panel) where the rms amplitude (^/o^) is shown as a function of 
count rate. To produce this plot the individual rms estimates (Fig.|6| 
panel 3) were sorted by count rate and binned by flux (such that 
there were 20 estimates per bin). The error on the mean rms was 
calculated in the standard fashion (see above). This indicates that 
the source does show a form of genuine non-stationarity: the abso- 
lute rms amplitude of the variations increases, on average, as the 
source flux increases. This effect has been noted in other Seyferts 
(Uttley & M'^Hardy 2001; Edelson et al. 2002; Vaughan et al. 2003) 
and is due to a linear correlation between rms and flux (see Uttley 
et al. in prep, for further discussion of this effect). Non-stationarity 
of this form can be 'factored out' by using the normalised ampli- 
tude (Fvar or (T^xs) instead of the absolute values. Normalising 
each variance (or rms) estimate by its local flux removes this trend. 
The bottom panel of Fig.|9|shows that Fvar is indeed constant with 
flux (fitting a constant gave = 7.8 for 9 dof). Fig.|6|(panels 5 
and 6) shows JVar and its average as a function of time; the aver- 
age is consistent with staying constant (x^ — 5.8 for 9 dof). The 
variability of Mrk 766 does show genuine (strong) non-stationarity, 
in the sense that the absolute rms increases linearly with flux, but 
this trend can be removed by using normalised units - JVar (and 
therefore the normalised excess variance) is consistent with being 
constant (with time and flux). 

The above analysis suggests that, after accounting for the ef- 
fect of the rms-flux correlation, there is no other evidence for 
strong non-stationarity in the rapid variability of Mrk 766. This 
was confirmed using the s-test of Papadakis & Lawrence (1995; see 
their Appendix A). A periodogram was calculated for three consec- 
utive segments of 3.4 x 10* s duration, and normalised to fractional 
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Figure 9. Top panel: The average absolute mis amplitude (-^ ""xs^ ^ ^ 
function of flux for the Mrk 766 light curve (compare with Fig.|^. Bottom 
panel: The average it fractional rms amplitude (•\/ ""nxs^ ^ function 
of flux. Clearly the absolute rms amplitude is a function of flux, but this 
dependence is removed in the fractional rms. 



units (see Appendix |aJ. The s value was computed by comparing 
periodograms at frequencies below 2 x 10^'^ Hz (above which the 
Poisson noise power becomes comparable to the source variabil- 
ity). For each pair of periodograms the value of s was within the 
range expected for stationary data (specifically |s| < 1, within one 
standard deviation of the expected value). 



7.3 rms spectrum 

The variability amplitude as a function of energy was calculated 
by measuring Fvar from light curves extracted in various energy 
ranges. The results are shown in Fie. 1101 and the errors were cal- 
culated using equation IB 21 to account for the effect of Poisson 
noise. The variability amplitude is clearly a function of energy, 
i.e. Mrk 766 shows significant spectral variability. This was con- 
firmed by a Fourier analysis of the light curves in different energy 
bands (Vaughan & Fabian 2003) which revealed complex energy- 
dependent variability. 

The rms spectrum was re-calculated for light curves contain- 
ing only single pixel (pattern 0) events and again for double pixel 
(patterns 1^) events. These two sets of data were extracted from 
the same detector and using identical extraction regions etc. Af- 
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Figure 10. rms spectrum of Mrk 766 measured using EPIC pn light curves 
with 1000 s bins. 



ter accounting for the difference in count rate between single and 
double pixel events the two sets of light curves should be identi- 
cal except for the effects of Poisson noise. The two rms spectra 
should be the same except for Poisson errors. Comparing the ra- 
tio of the two rms spectra using a x'^-test (against the hypothesis 
of unity ratio) gave = 25.3/20 degrees of freedom. Compar- 
ing the difference of the two rms spectra to the hypothesis of zero 
difference gave identical results and shows the two rms spectra are 
indeed fairly consistent. This test indicates that for real data the er- 
ror formula given above does provide a reasonable description of 
the uncertainty induced by photon noise. 



8 DISCUSSION 

The analysis of stochastic processes, such as X-ray variability of 
AGN, is conceptually different from the analysis of deterministic 
data such as time-averaged spectra (see discussions in e.g. Jenkins 
& Watts 1968; Priestley 1981; Bendat & Piersol 1986; Bloomfield 
2000). For example, when observing the spectrum of a constant 
source one expects repeatability of the data to within the limits set 
by measurement errors, i.e. each new realisation of the spectrum 
should be consistent within the errors. In AGN variability analy- 
sis it is the signal itself that is randomly variable; one does not 
expected repeatability of quantities such as the mean or variance. 
These statistical moments will change (randomly) with each new 
light curve even if there are no measurement errors. 

The stochastic nature of red noise processes means that it is 
usually only their average properties that can provide physical in- 
sight. Non-deterministic data should be handled statistically. For 
example, it is customary to examine the timing properties of X-ray 
binaries using PSDs estimated from the average periodogram of an 
ensemble of light curves (e.g. van der Klis 1995). Averaging over 
many independent realisations reduces the random fluctuations in- 
herent in the noise process. 

In most AGN timing studies however there are rarely enough 
data to construct averages in this way (but see Papadakis & 
Lawrence 1993 and Uttley et al. 2002 for more on PSD estima- 
tion for AGN). As a result of this relative lack of data, AGN tim- 
ing studies often emphasise the properties of a single light curve. 
But emphasis on the detailed properties of any single realisation 
of a stochastic process can be misleading. For example, AGN light 



curves show large fluctuations in variance. These changes provide 
little insight as they are expected even when the underlying physical 
process responsible for the variability is constant. Rather, they may 
simply be statistical fluctuations intrinsic to the stochastic process. 
All red noise processes show random fluctuations in both mean 
and variance and the variance will be distributed in a non-Gaussian 
fashion with a large scatter (see sectionsU2|and|5}. 

Previous claims of non-stationary variability based on changes 
in variance (e.g. Nandra et al. 1997; Dewangan et al. 2002; Gliozzi, 
Sambruna & Eracleous 2003) should therefore be treated with cau- 
tion since they did not account for this intrinsic scatter (see also 
section 3.3.1 of Leighly 1999 for a discussion of this point). Real 
changes in the PSD would indicate genuine non-stationarity and re- 
flect real changes in the physical conditions of the variability pro- 
cess. Such changes can be measured from the average properties of 
the light curve, such as the averaged periodogram or the averaged 
variance (see section|5j. 

A different issue is that differences between the variance of si- 
multaneous light curves obtained in different energy bands can be 
examined using the excess variance (or Fvar) statistic. It is possible 
to estimate the uncertainty in the excess variance due to errors in 
the flux measurements. This uncertainty, accounting only for mea- 
surement (e.g. Poisson) errors, can be used when testing for spectral 
variability, as demonstrated in sectionQ 

Estimators such as the excess variance provide a useful, if 
crude, means of quantifying the variability of AGN. Even though 
the stochastic nature of AGN light curves makes it difficult to es- 
timate variability amplitudes robustly from short observations, the 
excess variance can provide useful information. For example an 
analysis of the excess variances measured from short observations 
of Seyfert 1 galaxies demonstrated that the variability amplitude (at 
a given range of timescales) is inversely correlated with the lumi- 
nosity of the source (Nandra et al. 1997; Leighly 1999; Markowitz 
& Edelson 2001). Although random fluctuations in variance are ex- 
pected for AGN light curves the range of variances observed is far 
larger than could be accounted for by this effect alone. Another 
example is given in section lT^ when it is demonstrated that the av- 
erage variance of Mrk 766 is a function of the flux of the source. 
A similar effect has been observed in X-ray binaries (Uttley & 
M'^Hardy 2001). A discussion of the implications of this result will 
be given in Uttley et al. in prep. 



9 CONCLUSIONS 

This paper discusses some aspects of quantifying the variability of 
AGN using simple statistics such as the variance. Various possible 
uses of these are presented and some possible problems with their 
significance and interpretation are brought to light. The primary 
issues are as follows: 

(i) In order to search for non-stationary variability in an ensem- 
ble of short light curves (or short light curve segments) one can test 
whether the individual variances are consistent with their mean. 
Two practical methods are presented (sections i5.2.1l and l5.2.2t . 

(ii) In the first method the mean variance and its error are cal- 
culated at various epochs by binning the individual variance esti- 
mates. This is most useful when searching for subtle changes in 
variability amplitude but requires large datasets (in order that the 
variance can be sufficiently averaged). 

(iii) In the second method the individual variance estimates are 
compared with the expected scatter around the mean. The expected 
scatter is calculated using Monte Carlo simulations of stationary 
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processes. The table gives some examples of the scatter expected 
for various PSD shapes typical of AGN. This table can therefore 
be used to provide a 'quick look' at whether the observed fluctu- 
ations in the variance are larger than expected. One drawback is 
that, because the intrinsic scatter in the variance is rather large for 
red noise data, this method is only sensitive to very large changes 
in the variability amplitude. Another drawback is that one has to 
assume a shape for the PSD. 

(iv) The excess variance can also be used to quantify how the 
variance changes as a function of energy (section l6^ . An approx- 
imate formula is presented (based on the results of Monte Carlo 
simulations) that gives the expected error in the excess variance 
resulting from only observation uncertainties (flux errors such as 
Poisson noise). This can be used to test for significant differences 
in variance between energy bands. If the normalised excess vari- 
ances (or FvarS) are found to differ significantly between energy 
bands this implies the PSD is energy dependent and/or there are 
independently varying spectral components. 

(v) Possibly the most robust yet practical approach to variabil- 
ity analysis from AGN data is to test the validity of hypotheses 
using Monte Carlo simulations. This approach has yielded reliable 
PSD estimates for Seyfert galaxies (Green et al. 1999; Uttley et al. 
2002; Vaughan et al. 2003; Markowitz et al. 2003) and has been 
used to test the reliability of cross-correlation results (e.g. Welsh 
1999) amongst other things. Section|3|discusses some methods for 
simulating red noise data. 
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APPENDIX A: PERIODOGRAM NORMALISATION 

The periodogram is calculated by normalising the modulus-squared 
of the DPT (see equation 0: 

P{f,)^A\DFT{f,)\\ (Al) 

There are a variety of options for the normalisation A used 
in the literature, each has desirable properties. In the mathematical 
literature on time series analysis a normalisation of the form A = 
2/N is standard (e.g. Priestley 1981; Bloomfield 2000). However, 
this normalisation is generally not used for time series analysis in 
astronomy because the periodogram then depends on flux of the 
source and the binning of the time series. Below are listed three 
of the most commonly used normalisations, which only differ by 
factors of x, the mean count rate in cts/s (Aabs = xA-L^aj^y = 
x'^ A-^-^^i ). The factor of two is present in all these normalisations 
to make the periodogram 'one sided,' meaning that integrating over 
positive frequencies only yields the correct variance. 
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(i) ^,^,2 = 2ATsamp/S^A^ - defined by van der Klis (1997) 
(see also Miyamoto et al. 1991). This is the normalisation most 
often used in analysis of AGN and X-ray binaries because the inte- 
grated periodogram yields the fractional variance of the data. The 
units for the periodogram ordinate are (rms/mean)'^ Hz~^ (where 
rms/mean is the dimensionless quantity -FVar), or simply Hz~^. 

If a light curve consists of a binned photon counting signal (and 
in the absence of other effects such as detector dead-time) the ex- 
pected Poisson noise 'background' level in its periodogram is given 
by 

_ 2[X + B) Ar,an.p ,,,, 

where x is the mean source count rate, B is the mean background 
count rate, ATsamp is the sampling interval and ATbin is the time 
bin width. The factor of ATsamp /ATbin accounts for aliasing of 
the Poisson noise level if the original photon counting signal con- 
tained gaps. If the light curve is a series of contiguous time bins 
(i.e. ATbin = ATsamp) and has zero background (which is approx- 
imately true for many XMM-Newton light curves of AGN) then this 
reduces to Pnoisc = 2/x. 

For a light curve with Gaussian errors (Jcrr,* the noise level in the 
periodogram is 



(A3) 



ATbin 



(ii) ^Lcahy = 2 ATsamp /^A'^ - Originally due to Leahy et al. 
(1983). This has the property that the expected Poisson noise level 
is simply 2 (for continuous, binned photon counting data). If the 
light curve consists only of Poisson fluctuations then the peri- 
odogram should be distributed exactly as X2 ■ It is this property that 
makes this normalisation the standard for searching for periodic 
signals in the presence of Poisson noise (see Leahy et al. 1983). 
If the input light curve is in units of ct s~^ then the periodogram 
ordinate is in units of ct s~^ Hz~^. 

(iii) Aabs ~ 2 ATsamp /A' - this is the normalisation used in 
equation |2| This gives the periodogram in absolute units [e.g. (ct 
s~^)^ Hz~^] and so the integrated periodogram gives the total vari- 
ance in absolute units [e.g. (ct s~^)^] For a contiguously binned 
light curve with Poisson errors the noise level is Pnoisc = 2a;, and 
for Gaussian errors the noise level is Pnoisc = 2ATbincr|rr. 



APPENDIX B: MONTE CARLO DEMONSTRATION OF 
POISSON NOISE INDUCED UNCERTAINTY ON EXCESS 
VARIANCE 

To estimate the effect on a^xs due only to Poisson noise the basic 
strategy was as follows. 

(i) Generate a random red noise light curve. This acts as the 
'true' light curve of the source. 

(ii) Add Poisson noise, i.e. draw fluxes from the light curve ac- 
cording to the Poisson distribution. This simulates 'observing' the 
true light curve. Error bars were assigned based on the 'observed' 
counts in each bin (V counts). 

(iii) Measure the normalised excess variance (Tnxs °f ob- 
served light curve. This will be different from the variance of the 
true light curve because of the Poisson noise. 

Steps 2 and 3 were repeated, using the same true light curve, to 
obtain the distribution of cr^xs- Fie. lBll shows some results. In this 
example the 'true' light curve was generated with a PSD and 



normalised to a pre-defined mean and variance, e.g. S'^/x^ = 0.04 
(Pvar ~ 20%). This light curve was then observed (i.e. steps 2 
and 3 were repeated) 10* times**. The three panels correspond to 
different mean count rates for the true light curve (i.e. different S/N 
of the observation). The (la) widths of the cr^xs distributions are 
Monte-Carlo estimates of the size of the error bars on cr^xs due to 
Poisson noise. 

As is clear from Fig. IBll the distribution of (Tnxs becomes 
narrower, i.e. the error on cr^xs g^ts smaller, as the S/N of the 
data increases. Obviously in the limit of very high S/N data the 
measured value of a^xs t^nd to the 'true' value (in this case 
0.04), i.e. err((T^xs) — > as counts — > 00. It should also be 
noted that the distributions are quite symmetrically centred on the 
correct value, indicating that (Tnxs an unbiased estimator of the 
intrinsic variance in the light curve, even in relatively low S/N data. 

In order to assess how the error on a^xs changes with S/N, 
the width of its distribution was measured from simulated data at 
various different settings of S/N and intrinsic variance (i.e. S'^ /x"^). 
Width of the distribution at each setting was calculated from only 
500 'observations' of each light curve. In order that no particular 
realisation adversely affect the outcome, and to increase the statis- 
tics, this was repeated for 20 different random light curves (of the 
same fractional variance) and the width of the cr^xs distributions 
were averaged (i.e. the whole cycle of steps 1-3 was repeated 20 
times). Thus for each specified value of S/N and fractional variance, 
the error on cr^xs is estimated from 10* simulated 'observations.' 
These Monte Carlo estimated errors on the normalised excess vari- 
ance are shown in Fie. lB2l 

The solid lines show the functions defined by equation 1111 
(which was obtained by fitting various trial functions to the Monte 
Carlo results). Clearly this equation gives a very good match to the 
Monte Carlo results. 

If the variability is not well detected, either because the S/N 



is low or the intrinsic amplitude is weak, then 5* 



It is the 



first term on the right hand side of equation fTTI that dominates. If 
the variability is well detected, i.e. ^ (Terr, then it is the second 
term that dominates. 



err((jNxs) 




O ~ Ocrr 



(Bl) 



In the former case the deviations from the mean are dominated 
by the errors and the fluxes are approximately normally distributed. 
In this regime the error equation becomes the same as that given 
in equation A9 of Edelson et al. (2002). In the latter case the de- 
viations in the light curve are enhanced by the intrinsic variance. 
The second term is similar to the first except multiplied by a factor 



2crxs/cr|rr to account for this. 



err(P, 



Eguation ll ll can be used to give the uncertainty on F^, 
1 



2Pv: 



-err(CTNXs) 



\ 




thusly 



(B2) 



As this measures only the effect due to Poisson noise, the results are 
largely independent of the details of the Hght curve, including the PSD, as 
long as the flux is non-zero throughout the Hght curve. This was confirmed 
by repeating the above experiment using data produced from PSD slopes in 
the range q = — 2. 
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Fi gure B2. Width of the distribtition of (^^xs (i^^^tilting from Poisson 
noise) as a function of the number of counts per bin. Compare with Fig. lBll 
The solid curve shows the function described in the text (eauation ll II . 



errors on the fluxes. It does not account for the intrinsic scatter in 
the fluxes inherent in any red noise process. 

This paper has been typeset from a TgX/ file prepared by the 
author. 




of 

the same light curve. In each case the 'true' f^xs ^'^'^ (dotted line). The 
top panel used the lowest S/N data, the bottom panel used the highest S/N 
data. The mean number of counts per bin in the simulated light curves was 
15 (top), 30 (middle) and 100 (bottom). As the S/N increases (count rate 
increases) the distribution of o"^xs becomes narrower. (Note: this is differ- 
ent from Fig.|3] which shows how the variance changes between different 
realisations of the same stochastic process.) 



and this is the equation used to derive the errors shown in Fig. llOl 
In the two regimes this becomes: 




(B3) 



In the first instance, when the variability is not well detected, cr^xs 
should be preferred over Fvar as negative values of cr^xs possi- 
ble. Additional Monte Carlo simulations confirmed the above equa- 
tions are valid for both Gaussian and Poisson distributed flux errors. 
It is worth reiterating that this error accounts only for measurement 



